Absorbing processes in Richardson diffusion: analytical results 
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We consider the recently addressed problem of a passive particle (a predator), being the center 
of a "sphere of interception" of radius R and able to absorb other passive particles (the preys) 
entering into the sphere. Assuming that all the particles are advected by a turbulent flow and that, 
in particular, the Richardson equation properly describes the relative dispersion, we calculate an 
analytical expression for the flux into the sphere as a function of time, assuming an initial constant 
density of preys outside the sphere. In the same framework, we show that the distribution of times 
of first passage into the sphere has a t~ 5 ^ 2 power law tail, seen in contrast to the t~ 3 ^ 2 appearing in 
standard 3-D diffusion. We also discuss the correction due to the integral length scale on the results 
in the stationary case. 
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The statistics of relative dispersion of scalars in a turbulent flow is a topic of great theoretical and practical interest. 
It is well known, since the pioneering work of Richardson in 1926 1 , that relative particle dispersion in a turbulent 
i— !■ flow is superdiffusive, in particular that the mean-square separation between the particles grows with the third power 
j of the time. Richardson's explanation for this phenomenon was that, as the distance between particles increases, the 
effective diffusion should also increase since larger and larger eddies are involved. He was able to measure that the 
effective diffusion is proportional to r 4 / 3 , where r is the particle separation, and suggested a 3-D diffusion equation 
for the interparticle distance (now known as the Richardson equation) with a diffusion coefficient growing like r 4 / 3 
that indeed reproduces the correct superdiffusive behavior. Apart from practical implication, Richardson's result may 
i " 1 \ be considered a posteriori as the background for Kolmogorov's development of the idea of universality in turbulent 
flows^. 

Indeed, within Kolmogorov theory, the 4/3 exponent is simply obtainable by a dimensional argument: we will 
assume in the following that both the initial and the final separation are much smaller than the integral length scale 
L (the typical scale of the largest eddies) and much larger than the Kolmogorov scale r\ = (i^/e) 1 / 4 , where v is the 
viscosity and e is the mean dissipation rate of turbulent kinetic energy. Then, the effective diffusion can depend only on 
' the interparticle distance r and the strength of the field, which is proportional to r 1 / 3 . The result is -D e // = Ce 1 / 3 ?* 4 / 3 , 
where C is a universal constant. Plugging this coefficient into a 3D diffusion equation in spherical coordinates, one 
\& obtains: 
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d t f(r,t) = V P [D e //(r)V r /(r,t)] = Ce 1 / 3 ^ 2 ^ 10 / 3 ^/^, t)]. (1) 
By solving the equation with the intial condition r(t — 0) = 0, one can easily show that: 
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> (r 2 (*))cxei 3 . (2) 

Nowadays, we have a far better comprehension of turbulent phenomena, and we know that Kolmogorov theory is only 
approximatively valid because of intermittency. Still, the predictions of Eq.JSJ are found to be practically unaffected 
by intermittency corrections^: the robustness of this result is related to the fact that the mean dissipation rate e 
appears linearly into the equation. 

Actually, many phenomena of interest in chemical and biological sciences involving diffusion fall in the class of 
the absorption problems. A typical example are the so-called reaction rate problems in chemistry- : consider a 
particle A surrounded by an uniform distribution of diffusing particles B 7 and able to react instantaneously with 
them, A + B — > A, when they fall inside a sphere of radius R, centered in the position of A. The problem is 
then to calculate the flux inside the sphere as a function of time. In 3D Brownian diffusion, the time dependent 
flux writes <&(i) = AttRV[1 + R/y/nTH] where V is the diffusion coefficient; this classical result was firstly derived 
by Smoluchowski 6 . The same problem in the 2D case, formulated as the mean area spanned by a disc undergoing 
Brownian motion in a time t, was solved by Kolmogorov and Leontovitch some years lateri Another problem one 
could consider is the distribution of the first passage times: if there are few reactants, then it may be interesting to 
know the probability for a reaction to occur in a time t. 

The reaction rate problem was recently addressed for the description of plankton predator-prey dynamics in turbu- 
lent flows&S. Here the particle A (the predator) belongs to a plankton species feeding on the preys B, and one would 
like to calculate the amount of preys consumed per unit time. The biological assumptions underlying this model are 
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essentially that the predator is able to eat instantaneously every prey in its "sphere of interception" of radius R, and 
that both the predator and the preys are passive particles, i.e. their velocities are exactly the velocity of the fluid at 
their positions. 

On the physical side, it is tempting to face these absorption problems in turbulent fluids using a simple approach 
based on the Richardson equation; still, it is not trivial to assess whether this effective description of the evolution 
of the interparticle distance is valid when one considers problems that are more subtle than growth of the the mean 
square separation with time. We just point out here that some of these problems (like the first passage times problem) 
are known to be related to backward diffusion, and recentljiiS it has been shown that non-Gaussian tails and temporal 
correlations in the Eulerian flow may imply relevant differences between the forward and backward diffusion, like a 
different, non-universal value of the constant C. Still, it is encouraging that some predictions based on the Richardson 
equation for absorbing processes have been tested experimentallySiAi and numerically 12 . 

In this Letter, we show how these problems can be analytically tackled assuming the validity of the Richardson 
equation. Our main results are the analytical expression of the time-dependent flux, that is the solution of the reaction 
rate problem for the Richardson diffusion, and the calculation of the exponent of the power-law tails of the first passage 
time distribution. We also discuss, in the stationary case, the corrections due to the presence of an integral length 
scale L above which the diffusion coefficient does not depend on r anymore. 

Eq. is casted into a standard Fokker-Planck formic by defining the radial density p(r,t) — r 2 f(r,t): 



d t p(r,t) = Ce^ 3 d r 
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r 1 / 3 ^, t) + d r r 4 ^p{r, t) 



In order to simplify the expression, we make the substitution r 1 / 3 
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d t g(x,t)^Ce^ 3 d x 



g{x,t) + d x g{x,t) 



(3) 



(4) 



Notice that this equation is effectively the same as the equation for a 9-dimensional diffusion in radial coordinate. 
The subject of our study is obtained by taking the adjoint of the operator in the r.h.s of Eq.QJ, that is the backward 
Kolmogorov equation corresponding to the Richardson equation: 



d t f(x,t) = Ce x ' 3 



-d x f( x ,t) + d 2 j(x,t) 



(5) 



Eq.jnj) constitutes the appropriate equation to describe absorbing processes; indeed, by imposing absorbing bound- 
ary conditions on Eq. JSJ, f(R 1/>3 ,t) = /(oo,i) = Vi, the function f(x,t) represents the probability of not leaving 
the interval [i? 1 / 3 ,^) in a time t, i.e. of not being absorbed into a sphere of radius R, for a particle starting at 
position x. It can also be shownM that, calling 8(x, t) the radial concentration of reacting chemicals (or preys, in the 
language of reference^), and imposing slightly modified boundary conditions 9(R 1 ^ 3 ,t) = 0, #(oo,i) = p where p is 
the initial (constant) density of preys, then the function 8(x,t) still evolves according to Eq.©. 

By introducing the Laplace transform of f(x,t) with respect to time, f(x,t) — dte~ st f(x,t), Eq. JSJ) becomes: 



-d x f(x, s) + d 2 J(x, s) - sf(x, s) = 0. 



(G) 



If we find a solution h(x) of Eq. JH| for s = 1 which is real and positive for x > 1, then it is easy to show that 

s /i(Vs) 



(7) 



is the complete solution to Eq. I© with the boundary conditions corresponding to the absorption problem, while the 
reaction rate problem is solved by 9(x, s) = p[l/s — f(x, s)]. 

The function h(x) is, up to a multiplicative constant, found to be 



h{x) = x 7/2 K 7/2 (x) = y|^ 7 (15 



15x + 6x 2 + x 3 )e~ x 



(8) 



where K is the modified Bessel function of the second kind^, and we have written its expression in terms of simple 
functions. 
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Let us now move to the calculation of the time-dependent flux into the sphere: 



<PM*) = --D(r) — g- — \r=R 



and the integrated flux $«(s) = AnR 2 (j)ji(s) equals 
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s K 7/2 (B}/3y/s/(CeW)) 
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= 47r ^ 
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Next we apply the inverted Laplace transform. Notice that the first term in Eq. I jlOj l corresponds to the t 
steady state flux: = 28nCe 1 ^ 3 R 7 ^ 3 /3. We invert the Laplace transform of the following function: 



1 K 5/2 (^) 



3 + 3^ + 3 



v^K 7/2 (y/s~) 15 + 15\/i+12s + s 3 / 2 
and the inverse Laplace transform is given by the Bromwich integral 
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ds /(s) e ts ds 



(10) 
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where a is an arbitrary real number such that all the singularities of the integrated function lie on the left of the 
contour of integration. To avoid problems with branch cuts, we make the substitution ^fs = z 
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1 r 6z + 6z 2 + 2z 3 

2t7i J Tl Z 15+15z + 12z 2 + 



dz 



(13) 



where the contour path is Ti in the picture. We evaluate the integral over the closed contour of the figure: it is 
possible to calculate the poles of the integrand, they have negative real parts and are thus outside the integration 
contour. It is also possible to show that the contribution of the arcs vanishes according to Jordan's Lemma. Then we 




FIG. 1: Integration contour in the complex plane. 



have: 
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dz 2zf(z 2 )e tz dz 
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dz 2zf(z 2 )e t 
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~2ir]_^ dy 225 - 135t/ 2 + 114y 4 + t/ 6 ' 



e v dy 



(14) 
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where the variable is y — iz. The ratio of polynomials in l|14fl is obtained as the real part of the function 2zf(z 2 ) 
calculated on the imaginary axis and its imaginary part vanishes when integrated since it is an odd function of the 
variable y. 

Going back to our problem, the time solution is found to be: 



28irCe 1 / 3 pR 7 / 3 



2Ce 1 l 3 pR 1 ' 3 



dy 



2y 4 (i8 + y 2 ) 



225 - 135y 2 + lUy 4 



CeV 3 « „3 
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and, as it is expected, can be cast into the scaling form &R(t) = Cpe 1 / 3 R 7 / 3 ui(Cc 1 / 3 t/R 2 / 3 ). The shape of the scaling 
function w(r), where the non-dimensional time r is equal to Ce 1 / 3 t/ R 2 / 3 , is shown in Fig.(J2|. Notice also that the 
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FIG. 2: Plot of the scaling function lj(t) 

flux diverges like r -1 / 2 when r — > 0. 

We move now to the problem of estimating the probability for a particle to be absorbed in a time t, given the 
starting distance x. We already know the Laplace transform of this probability f(x,s), which is the solution of 
Eq.©. Our goal is to calculate the inverse transform: 



f(x,t) 
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f(x, s)e ts ds 
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In order to do that, we call z = y/s and perform the integral over the same contour as in the previous problem. In this 
case, we have to avoid a singularity in the origin which gives a time-independent contribution corresponding to the 
finite probability of being absorbed for t — > oo. The time behavior is then given by the principal part of the following 
integral: 
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k(x,y) 
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with: 



k(x, y) = 6y(75(ir - 1)) - 5(-l + Zx{2 + x(5x - 2)))y 2 + 2.t 2 (15o: - l)y 4 cos(yx - y) + 
+6(75 - 15(.t - 2) (2a; - l)y 2 - x(5 + 3x(25x - A))y A + 5x 3 y 6 ) sin(y - xy) . 



(18) 



The first non-zero term in the expansion is proportional to t~ 3 / 2 , meaning that the first-passage time have a t~ 5 / 2 
power law tail. In Fig. we compare simulations of the first passage process with this theoretical prediction, by a 
simple numerical integration of the Langevin equation corresponding to the Fokker Planck equation J3J| . 

Now, we want to show that the effect of the integral length scale is just to introduce a small correction to the 
stationary flux, by solving the stationary problem V,.[Z3 e //(r)V r /(r)] = 0, where D e ff(R) is equal to Cer 4 / 3 for 
r < L, being L is the integral length scale defined in the beginning, and it is constant, D e tt(r) — CeR 4 / 3 , for 
r > L. In 3D, the equation becomes r~ 2 D e ff(r)d r f(r) — const, where the value of the constant is determined by the 
boundary conditions /(oo) = p and f(R) = 0: 



const = p 



dr 



2 Deff(R) 
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Ce 1 / 3 (3i?- 7 /3 + 4L-V3) 



(19) 
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FIG. 3: Probability density function of the first passage times. The statistics for different initial conditions, indicated in the 
legend, is collected from 5 • 10 6 realization of the Langevin equation corresponding to eq. A t~ 5 ^ 2 power law is also shown 
for reference. 



so that the stationary flux becomes: 

^ , . 28 7 rCe 1 /3 pi? 7/3 

= 3 + 4(W/3 ( 2 °) 

As one could expect, the effect of the integral length scale is to decrease the flux, and the correction can be safely 
neglected for realistic values of the parameters (R is supposed to be much smaller than L) . 

To conclude, by assuming the validity of a description a la Richardson for backward relative diffusion in turbulent 
flows, we show how simple problems like the calculation of the flux into a moving absorbing sphere and the first 
passage times distribution are analytically solvable. We also discuss the effect of the integral length scale on these 
results. Our solution could provide a basis for building up more realistic and detailed models of absorption phenomena 
in turbulent flows. 
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